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Abstract. We present Monte Carlo simulations of the dynamical evolution of 

,.1^ , the Oort cloud over the age of the Solar System, using an initial sample of one 

Mh' million test comets without any cloning. Our model includes perturbations due 

to the Galactic tide (radial and vertical) and passing stars. We present the first 

5^ . detailed analysis of the injection mechanism into observable orbits by comparing 

"ti ' the complete model with separate models for tidal and stellar perturbations alone. 

C^ , We find that a fundamental role for injecting comets from the region outside the 

loss cone (perihelion distance g > 15 AU) into observable orbits (q < 5 AU) is 

played by stellar perturbations. These act in synergy with the tide such that the 

total injection rate is significantly larger than the sum of the two separate rates. 

^^ • This synergy is as important during comet showers as during quiescent periods and 

VQ ' concerns comets with both small and large semi-major axes. We propose different 

1/^ , dynamical mechanisms to explain the synergies in the inner and outer parts of the 

^S| ' Oort Cloud. We find that the filling of the observable part of the loss cone under 

l ' I normal conditions in the present-day Solar System rises from < 1% for a < 20 000 AU 

^ ! to about 100% for a ^ 100 000 AU. 

^— V ^ Keywords: Galactic dynamics, Monte Carlo simulation, long period comets, Oort 

cloud 
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1. Introduction 

When analyzing the distribution of original inverse semi-major axes 
of long-period comets, Oort (1950) concluded that the near-parabolic 
spike of this distribution reveals a distant reservoir of comets (the 
'Oort Cloud'). His favoured mechanism of injection of comets from this 
reservoir into observable orbits {i.e., with small perihelion distances) 
was the passage of stars in the vicinity of the reservoir, whereby the 
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long-term reshuffling of angular momenta would ensure a continuous 
infeed into the innermost part of the Solar System. 

Until the 1980's stellar perturbation was the only mechanism con- 
sidered, when issues concerning the injection of comets from the Oort 
Cloud were discussed {e.g., Rickman, 1976; Weissman, 1979; Fernandez, 
1980; Hills, 1981; Remy and Mignard, 1985). However, by that time 
it became clear that the tidal action of the Galaxy as a whole must 
also have an important influence - especially the part corresponding to 
the z-dependent disk potential (Heisler and TVemaine, 1986). This was 
verified by noting that the Galactic latitudes of perihelia of new Oort 
Cloud comets have a double-peaked distribution that is characteristic 
of the disk tide (Delsemme, 1987). 

An important paper by Duncan et al. (1987) treated the formation 
of the Oort Cloud and showed that the characteristic time scale for 
changing the perihelion distances, independent of the semi-major axis, 
is shorter for the Galactic disk tide than for the stellar perturbations. 
This has been verified by later analytic work, e.g., by Fernandez (2005), 
and further numerical simulations of Oort Cloud evolution {e.g., Heisler 
et al., 1987) have given support to the dominance of Galactic tides for 
comet injection. 

Consequently, stellar perturbations have come to be practically ne- 
glected as a source of comet injection - except when discussing "comet 
showers" (Hills, 1981) arising from rare stellar passages through the 
dense, inner parts of the Oort Cloud. The importance of stellar per- 
turbations for randomizing the orbital distribution of the Oort Cloud 
and thus keeping the relevant infeed trajectories of the disk tide pop- 
ulated over long time scales has been realized (see Dybczyhski, 2002 
and references therein), but the actual injection is often seen as due 
only to the tide. Hence it should be limited to semi-major axes large 
enough for the tidal perturbation to bring the cometary perihelion at 
once from outside the "Jupiter-Saturn barrier" {i.e., perihelion distance 
q ^ 15 AU) into the observable region {q < 5 AU). The result is that 
one expects new comets to have aori ^ 3 x 10^ AU (Bailey and Stagg, 
1988; Levison et al., 2001; Fernandez, 2005). 

On the other hand, some recent papers indicate that this picture may 
have problems. The fractional population of the observable region - if 
fed only by Galactic tides - is small enough, and the orbital periods long 
enough, that the estimated total population of the Oort Cloud may be 
uncomfortably large (Charnoz and Morbidelli, 2007). And in addition, 
when non-gravitational effects are included into orbit determinations 
for new Oort Cloud comets (Krolikowska, 2006), the resulting original 
orbits tend to be of shorter periods, having smaller semi-major axes - 
often much smaller than 3 x 10^ AU. 
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Meanwhile, we have developed fast and accurate methods to treat 
both the Galactic tides (Breiter et al., 2007; Fouchard et al., 2007) and 
stellar perturbations (Rickman et al., 2005) in Monte Carlo simulations 
of Oort Cloud dynamics. This has allowed us to perform calculations, 
to be presented here, where the cloud is represented by as many as 
10^ sample comets and integrated over a time exceeding the age of 
the Solar System. This amounts to 5 x 10^^ comet-years of individual 
evolutions (or only slightly smaller due to the loss of comets during the 
simulation), which is much more than in all previous long-term Oort 
cloud simulations - e.g., 3 x 10^^ comet-years for Duncan et al. (1987), 
4 X 10^^ comet-years for Mazeeva (2006), and ~ 4 x 10^^ comet-years for 
Emel'yanenko et al. (2007) who used cloning. Heisler (1990) simulated 
~ 7 x 10^^ comet-years but only thanks to extensive cloning during the 
course of the simulation. In fact her long-term simulations (4.5 Gyr) 
concerned only ~ 10^ "tokens", i.e., comets actually followed, while 
these were meant to represent ~ 150 times as many comets by means 
of cloning. 

Our work is the first to study the mechanism of injection of comets 
from the Oort Cloud over the age of the Solar System by simulating 
and comparing different dynamical models. The reason why models 
involving both the Galactic tide and stellar perturbations gave a much 
higher flux of injected comets than those involving only stars (Heisler 
et al., 1987; Heisler, 1990) was never clarified, since comparisons with 
models involving only the tide were not made. In the present paper we 
concentrate on a comprehensive comparison of combined and separate 
models, thus describing and analyzing for the first time the synergy 
effect of Galactic tide and stars. 

We also take special care to define correctly the encounter velocities 
in our sample of passing stars, thereby arriving at somewhat larger 
values than those used previously. Finally, we study the relative filling 
of the observable part of the loss cone and the distribution of inverse 
semi-major axes of the injected comets. These studies are, however, 
only preliminary, since our current simulations do not include planetary 
perturbations, and thus we cannot account for those comets that arrive 
into the observable region after having "diffused" across the Jupiter- 
Saturn barrier in several revolutions. 

Our calculations are presented in Sect. 2, and in Sects. 3-5 we de- 
scribe our results in terms of the distribution of injection times into the 
inner planetary system, the flux of new, observable comets as a function 
of time, and the distributions of inverse semi-major axis and Galactic 
latitude of perihelion as well as loss cone filling at representative epochs. 
In Sect. 6 we discuss the results and summarize our conclusions. 
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2. Calculations 

As a simplifying assumption, we consider the Oort Cloud to have 
been formed instantaneously at a given epoch, and that its orbital 
distribution was isotropic to begin with. Thus the initial conditions are 
chosen with flat distributions of cos io, ujo, 0,o and Mq (we use common 
notations for the orbital elements, and the angles may be defined with 
respect to an arbitrary frame of reference). We consider a thermalized 
initial state of the cloud, where the semi-major axes {ao) are chosen 
in the interval 3 x 10^ < Oo < 1 x 10^ AU with a probability density 
oc a~^'^ (Duncan et al., 1987). The initial eccentricities (cq) are chosen 
with a density function oc Cq in such a way that the perihelia are outside 
the planetary system (g > 32 AU). We thus initialise 1 x 10^ comets. 

The Galactic parameters used for calculating the tidal effects are 
the same as described in earlier papers (Fouchard et al., 2007). The 
most important one for comparison with other investigations is the 
mid-plane disk density, which we take as 0.1 Mqpc~^ (Holmberg and 
Flynn, 2000). This is in agreement with Emel'yanenko et al. (2007), 
while Heisler (1990) used 0.18 Mqpc'^ (Bahcall, 1984). 

The simulation runs with a predefined set of 197906 stellar encoun- 
ters, occurring at random times during a lapse of tmax = 5 x 10^ yr, 
with random solar impact parameters up to dmax = 4 x 10^ AU, and 
with random stellar masses and velocities. Our procedure for creating 
each of these encounters is as follows. Let ^ denote a stochastic, random 
number with a uniform probability distribution on the interval [0, 1]. 

1/2 

The solar impact parameter is chosen as d = ^J x dmax, and the time 
of the encounter (specifically, the time of the star's perihelion passage) 
is t = ^i X tmax ■ The direction of stellar motion with respect to the Sun 
is defined in terms of Galactic latitude and longitude (6, t) such that 
sin 6 = 2^5 — 1 and i = ^ix 2tt, i.e., it has an isotropic distribution. The 
point at which the initial stellar velocity cuts the orthogonal impact 
plane is situated on a circle of radius d around the Sun, and its location 
is defined by an azimuthal angle (a) such that a = ^a >^ 27r. 

Next we choose the type of the star. We use 13 categories as in Rick- 
man et al. (2004) with parameters listed in Table I. To each category we 
associate one value of the stellar mass. These masses are generally taken 
as those of the archetypal spectral classes along the main sequence 
according to Allen (1985). However, in contrast to our earlier inves- 
tigations, we replace the archetypal mass of 18 Mq for BO stars by a 
weighted average of 9 Mq, considering that the less massive, later types 
(B2, B5) are much more common than the earlier ones. The relative 
encounter frequencies fi of Table I are taken from Garcia-Sanchez et al. 
(2001), where they were derived from the respective products of number 
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Table I. Stellar parameters. The types are mostly MK types for main 
sequence stars; 'wd' indicates white dwarfs, and 'gi' indicates giant 
stars. The encounter frequencies are given in number per Myr within 
1 pc. The following two columns list the solar apex velocity with re- 
spect to the corresponding type, and the spherical Maxwellian velocity 
dispersion. The last two columns give the mean heliocentric encounter 
velocity and its standard deviation according to our results. 



Type 


Mass 


Enc. freq. 


Vq 


O", 


(v-) 


c^v 




(Mo) 




(km/s) 


(km/s) 


(km/s) 


(km/s) 


BO 


9 


0.005 


18.6 


14.7 


24.6 


6.7 


AO 


3.2 


0.03 


17.1 


19.7 


27.5 


9.3 


A5 


2.1 


0.04 


13.7 


23.7 


29.3 


10.4 


FO 


1.7 


0.15 


17.1 


29.1 


36.5 


12.6 


F5 


1.3 


0.08 


17.1 


36.2 


43.6 


15.6 


GO 


1.1 


0.22 


26.4 


37.4 


49.8 


17.1 


G5 


0.93 


0.35 


23.9 


39.2 


49.6 


17.9 


KO 


0.78 


0.34 


19.8 


34.1 


42.6 


15.0 


K5 


0.69 


0.85 


25.0 


43.4 


54.3 


19.2 


MO 


0.47 


1.29 


17.3 


42.7 


50.0 


18.0 


M5 


0.21 


6.39 


23.3 


41.8 


51.8 


18.3 


wd 


0.9 


0.72 


38.3 


63.4 


80.2 


28.2 


gi 


4 


0.06 


21.0 


41.0 


49.7 


17.5 



density and mean velocity, ni{vi). A random number ^j is used to pick 
a stellar category i with the probability fi/J2fi- 

Finally, we choose the speed of the stellar motion in the following 
way. The velocity dispersions (o"*j) listed in Table I are taken from 
Garcia-Sanchez et al. (2001), and they correspond to the semi-axes 
of the velocity ellipsoids (aui, cr^i, a^i) listed by Mihalas and Binney 
(1981) using: cr^j = a^j + al^ + cr^j. For the pecuhar velocity {v^) 
of a star with respect to its LSR, we use a spherical Maxwellian by 
taking r/„, rj.^ and r/^ as three random numbers, each with a Gaussian 
probability distribution with expectance and variance 1, and com- 
puting v^ = a^i {(r/^ + 111 + Vw)/"^] • The star's heliocentric velocity 
is found by combining the vector v* with the Sun's peculiar velocity 
with respect to the star's LSR ("apex velocity") v©, whose absolute 
value is listed in Table I for each stellar category. We assume a random 
relative orientation of the two vectors and thus compute: 



V = {vli + vl- 2vQiV^ ■ C} 



1/2 



fl) 
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where C = cos 9 is taken as C = 2^c ~ 1; and 9 is the angle between 
the two vectors. 

Within each stellar category we have to account for the fact that 
the contribution to the encounter flux is proportional to V. Thus we 
define a constant, large velocity Voi = f0i + 3(T*i for each category, such 
that V is always smaller than Voj, and we take a new random number 
^v and keep the value just found for V, if ^v < V/Voi. Otherwise we 
repeat the computation of V until the ^y condition is fulfilled. This 
procedure was not followed in our previous investigations, leading to 
underestimates of the average stellar velocities. Further underestimates 
were caused by programming errors, and we caution the reader that the 
stellar velocities in Rickman et al. (2004; 2005) were systematically too 
small. This is clearly seen by comparing Fig. 1 of Rickman et al. (2005) 
with the data in Table I, which yield a mean velocity of 53 km/s with 
a dispersion of ~ 20 km/s. 

A few comments on the mean stellar encounter velocity are in order. 
The average peculiar velocity of stars in the solar neighbourhood is 
~ 40 km/s. This value was given by Hut and Tremaine (1985), and 
combining it in quadrature with a typical solar apex velocity of 23 km/s 
for the most common stellar categories (Table I), one gets a mean 
heliocentric velocity of ~ 46 km/s. Heisler et al. (1987) were the first 
to introduce the fiux-weighting into the selection of random velocities, 
but they neglected the solar apex velocity. In fact, their flux-weighting 
was somewhat different from ours, because they considered only one of 
the three velocity components, namely, the radial heliocentric velocity. 
But the encounter fiux is sensitive to the velocity with respect to the 
impact plane, i.e., the full speed of the star, instead of just the radial 
component. We have found that this difference has only a small effect 
on the resulting mean velocity, but we mention it for the sake of com- 
pleteness. In both cases, we find that the weighting raises the mean 
velocity by 6 — 7 km/s. This explains our mean velocity of 53 km/s as 
resulting from including both the solar motion and the flux-weighting 
procedure. Finally, let us compare with the mean encounter velocity of 
~ 46 km/s in Garcia-Sanchez et al. (2001). This resulted from a list of 
92 stellar encounters within 5 pc and 1 Myr of the present, compiled 
with the aid of Hipparcos data, but the authors showed that there was a 
serious bias against faint absolute magnitudes in this sample, affecting 
all stars with My > 4. Thus, the stars with the highest velocities 
were essentially lacking, and the resulting mean velocity is that of the 
inherently brighter, slower moving stars. 

Comparing with other investigators, we note that both Heisler (1990), 
Mazeeva (2006) and Emel'yanenko et al. (2007) based their stellar en- 
counter frequencies on the analysis by Heisler et al. (1987), who ignored 
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the solar motion - thereby underestimating the relative frequency of 
encounters with massive stars that have small velocity dispersions - and 
neglected the contribution of the massive giants. Our encounter sample 
contains as much as 3.5% of massive stars, i.e., the BO, AO, A5, FO, F5 
and 'gi' categories in Table I with an average mass of 2.3 Mq, while 
counting the stars in the absolute magnitude range that corresponds 
to this average mass in Heisler et al. (1987), one arrives at < 1% of 
the total encounter frequency. Since the massive stars have an average 
M/V ratio ^ 10 times larger than the red dwarf stars that dominate 
the sample, each such star will affect ~ 100 times as many comets. 
Hence one easily realizes that in our case a large fraction of the total 
stellar effect comes from the massive star category that is downplayed 
by the other investigators. This, to some extent, compensates for two 
other effects that make the stellar perturbations less efficient in our 
simulation. One is the larger encounter velocities, as already described, 
and the other is the total encounter frequency within 1 pc, which in 
our case is 10.5 per Myr, while for the others it is 13.1 per Myr. 

Our calculations of the heliocentric impulse imparted to the comet 
(at time t) are done using the Sequential Impulse Approximation, which 
guarantees a good accuracy at a low cost of computing time (Rickman 
et al., 2005). During the simulation we keep track of all the perihelion 
passages with their q values. At each perihelion a decision is taken about 
which method to use for the Galactic tide perturbation during the 
coming orbital period. The fastest method is a mapping (Breiter et al., 
2007) that analytically computes the orbital elements at the subsequent 
perihelion, but this is used only for elliptic orbits with semi-major axis 
less than a critical value that depends on the eccentricity (Fouchard 
et al., 2007). Otherwise we use numerical integration with a symplectic 
integrator for KS-regularized equations of motion (Laskar and Robutel, 
2001) in case 1/a > 10^^ AU~^ and the 15*^* order RADAU integrator 
(Everhart, 1985) for 1/a < 10"^ AU'^ 

During the orbital period in question, normally, several stellar pas- 
sages occur. On those occasions the osculating cometary orbit is subject 
to an instantaneous impulse. In the numerical integration regime for the 
Galactic tide, one always comes back to perihelion. But in the mapping 
regime, when the starting orbital period has elapsed, the comet may 
not be at perihelion because of the intervening stellar perturbations, 
and we then resort to numerical integration until the next perihelion 
passage takes place. 

The simulation proceeds for a maximum of 5 Gyr, unless an end state 
is reached. There are two such end states: either the comet reaches per- 
ihelion with q < qc= \f> AU (it is lost due to planetary perturbations). 
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or the comet reaches r = 4 x 10^ AU in outbound motion (it escapes 
directly into interstellar space). 

What we have just described is the full simulation of the "combined" 
model including both the Galactic tide and stellar perturbations. In 
addition, we have run two simulations that include only one or the 
other of the two dynamical effects. 



3. Injection time 

We will first consider the time needed to shift any cometary starting 
perihelion distance qo > 32 AU into a perihelion distance q < 15 AU. 
This is the time tinj required to inject a comet into the target zone, 
and we call it the injection time. We have thus scrutinized all three 
simulations, and for each injected comet in every simulation we noted 
its value of tinj. Let us now compare the statistics of injection times 
between all three dynamical models. 

The range of initial semi-major axes 3 000 < Oo < 100 000 AU is 
divided, according to logOo, into ten equal intervals. For each interval 
the following statistical parameters concerning the injection time are 
computed: its median value, its lower and upper quartiles (surpassed 
by 75%, respectively 25%, of the values), its lower and upper deciles 
(surpassed by 90%, respectively 10%, of the values), and finally its 
lower and upper percentiles (surpassed by 99%, respectively 1%, of the 
values). 

Figure 1 presents the comparison of tinj statistics by means of two 
plots. The left one (Fig. la) compares the model with only the Galactic 
tides to the one with only stellar perturbations, while in the right one 
(Fig. lb) the tides-only model is compared to the combined model. In 
each case we plot the statistical quantities versus Qq- The tides-only 
model is represented by filled squares, and for the other models we use 
half-filled circles. At any particular value of Oo, the symbols for each 
model are joined by vertical bars. A slight horizontal shift between 
the symbols has been introduced for easy distinction of the models, 
but the real Uq intervals are identical. The median values have been 
joined by curves (dotted for the tides-only model, and solid for the 
other ones). The grey dots show individual injections for the stars-only 
model (Fig. la) and the combined model (Fig. lb). 

For comparison with Duncan et al. (1987), who plotted a similar 
diagram (their Fig. 2), we include three lines. The one labelled Tg 
shows their "tidal torquing time" and the one labelled Tg shows the 
corresponding "stellar torquing time" , both as functions of a. These are 
meant to represent the typical time required to decrease the perihelion 
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Figure 1, Injection time versus initial semi-major axis. The semi-major axis range 
is split into ten equal intervals of log a. For each interval we plot different statistical 
parameters characterizing the distribution of injection times, as explained in the 
text. The dotted line labelled Tg corresponds to the tidal torquing time shown in 
Fig. 2 of Duncan et al. (1987), Ts to the corresponding stellar torquing time, and 
Tr to the period of {q,iG) oscillation imposed by the vertical tide. The left panel 
(Fig. la) compares the model with only tides to the one with only stars, and the 
grey dots are individual stellar injections. In the right panel (Fig. lb) the stars-only 
model is replaced by the combined model. 



distance from 25 AU to 15 AU in the two cases. The third line labelled 
Tr refers to the period of (g, ic) oscillation due to the tidal component 
normal to the Galactic plane (Duncan et al., 1987). 

By inspecting Fig. 1, we can make the following observations. First, 
compare the tides-only median curve with the T^ line. The two agree 
fairly well in the range from a ~ 6000 to 25 000 AU. This is natural, 
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because T^ is twice the average time it takes for the vertical tide to bring 
any Oort Cloud comet into the target zone, as long as it is on a relevant 
trajectory with qmin < 15 AU. For a ^ 6000 AU the median flattens out 
at about 2.5 Gyr, and this is due to the limit of our simulated interval at 
5 Gyr. Had we let the simulation run for a much longer time, we would 
have seen the median curve follow the T^ line to even smaller a values. 
For a ^ 25 000 AU we see how the median curve turns upwards, while 
Tr continues to decrease. This can be explained as a result of a quick 
stripping of comets from all trajectories with qmm < 15 AU, after which 
these have to be repopulated through the action of the non-integrable 
part of the tides. Since this works on a much longer time scale, it is 
obvious that the median of tinj has to increase. 

Already at this point we see evidence that the mean injection time 
- even in the tides-only model - does not follow the prediction of the 
Tr{a) dependence at all semi-major axes. Studying the Oort Cloud over 
a long enough time allows other parts of the Galactic tide than the 
simple, vertical component to take control of comet injections, at least 
in the outer parts of the cloud. But consider now the median curves of 
the two models that involve stellar perturbations. They are mutually 
quite similar, but they differ strongly from the tides-only curve except 
at a < 6000 AU. 

The mutual similarity - in spite of a much larger number of grey 
points (injections) in Fig. lb - means that the same basic mechanism 
is at work. We identify this as the angular momentum reshuffling by 
stellar perturbations. In Fig. la (stars-only model) this in itself makes 
comets diffuse all over angular momentum space so that some reach the 
target zone. In Fig. lb (combined model) the same angular momentum 
diffusion repopulates the "infeed trajectories" (with qmin < 15 AU) of 
the vertical Galactic tide, whereupon the comets are injected at a rate 
given by Tj. (a) . We interpret the flatness of the median curves at a level 
of roughly half the duration of our simulation as evidence that the time 
scale of angular momentum reshuffling is short enough to guarantee an 
injection rate that is at least as large during the first half (0 — 2.5 Gyr) 
as during the second half (2.5 — 5 Gyr). This is in agreement with 
the thermalization time scale reported by Duncan et al. (1987). The 
tendency for a slight decrease of the median ijnj at the largest semi- 
major axes is likely due to a progressive depletion of the outer parts 
of the cloud during the course of the simulation, while the reshuffling 
time scale is relatively short. 

We thus realize that the behaviour of the median injection time, 
generally speaking, has very little to do with any of the theoretical 
time scales. Let us now instead consider the lowest percentiles, since 
these give information on the quickest injections that - in principle - 
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might be governed by the Tg or Tg dependences. The lowest percentile 
of the tides-only model indeed decreases with Uq roughly parallel to 
the Tg line in the inner core of the cloud (a < 10 000 AU), but we 
nonetheless see a somewhat smaller slope. This tendency gets stronger 
with increasing Oq and finally turns into an increase outside 30 000 AU. 
We interpret this as due to the same repopulation of infeed trajectories 
by the non-integrable part of the tide as we discussed in relation to the 
median curve. 

The lowest percentile of the stars-only model shows a fall-off with ao 
that is interestingly slow in comparison with the Tg line. This appears 
to be related to the horizontal bands of grey points, which are cometary 
showers. For each ao interval in the inner core, the timing of the lowest 
percentile is that of the first shower reaching into that interval. The 
larger Oo, the sooner such a shower appears. But the showers also 
get weaker, being caused by more and more distant stellar encounters. 
Thus, in the outer parts of the cloud they are no longer of significance 
for defining the lowest percentile. Since this is instead controlled by 
a growing number of usual, inefficient stars passing through the outer 
regions, one has to wait longer. 

When we look at the lowest percentile of the combined model, we 
see that it follows the same decrease as the tides-only model in the 
inner core. Indeed, with a much larger number of injections, the comet 
showers have lost their importance, and as we shall see in Sect. 4, 
during the first Gyr the injections are largely controlled by the Galactic 
tides. But outside the inner core the lowest percentile now behaves with 
respect to that of the tides-only model in a similar way as the median 
does, and the reason is the same. Going to larger semi-major axes, in 
both models we see an increasing number of late injections, although for 
different reasons, and these determine the behaviour of most statistical 
parameters, causing them to decrease less rapidly than the Tg line. 

In summary we can state that we have found the theoretical time 
scales of Duncan et al. (1987) to give a rather poor representation of our 
statistics of injection times. On the other hand, we find evidence in the 
combined model for an important role being played by the repopulation 
of tidal infeed trajectories via stellar encounters - something that may 
be described as a synergy effect. This being said, one nonetheless real- 
izes that Tr is one of the basic time scales that govern this synergetic 
injection, the other one being the angular momentum reshuffling time 
scale of stellar perturbations. Let us now move to a discussion of the 
rate of injections and how this depends on time. 
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4. Time dependent injection flux 



The upper part of Fig. 2 shows a histogram plot of the number of 
comets injected into the observable region as a function of time from 
the beginning till the end of the simulation. Three histograms are shown 
together: the one in black corresponds to a model with only Galactic 
tides, and the grey one to a model including only stellar perturba- 
tions. Finally, the top, white histogram is for the combined model that 
includes both tides and stars. 
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Figure 2. The upper diagram shows the number of comets entering the observable 
zone per 50 Myr versus time. The white histogram corresponds to the combined 
model, the black histogram to the Galactic tide alone, and the grey histogram to 
the passing stars alone. The asterisks indicate the number of comets remaining in 
our simulation for the combined model at every 500 Myr with scale bars to the right. 
The middle diagram shows the excess number of injections into the observable region 
per 50 Myr in the combined model with respect to the sum of the stars-only and 
tides-only models. The lower diagram shows this excess expressed in percent of the 
mentioned sum. 

The first thing to note is the gradual decline of the injection flux 
of the tides-only model over a few Gyr, after which it stays at a very 
low level. The reason is clear. In the beginning, the phase space tra- 
jectories that in the regular dynamics imposed by the vertical tide will 
periodically lead into the "loss cone" [q < 15 AU) are populated just as 
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densely as any phase space domain and thus furnish an important flux 
of injections during the first period of (e, ic) oscillations. This amounts 
to a typical time scale of < 1 Gyr for much of the initial cloud, but the 
population on these trajectories is depleted by each entry into the loss 
cone, and there is no efficient way to replenish them without including 
stellar passages. On the longer time scale, we see only the feeble flux 
coming from (1) the infeed into the tidal injection trajectories by the 
non-integrable part of the tide; (2) the inner parts of the cloud, where 
the period of oscillation is very long (Fouchard et al., 2006). 

The other two histograms include the effects of stellar passages, 
and the stars are the same in both simulations. Therefore, we see the 
same comet showers appearing and the same quasi-quiescent periods 
in between. The white area at the top of each bin corresponds to the 
extra contribution of the combined model as compared with that of 
the stars only. If the numbers plotted in the white, grey and black 
histograms are called Nc-, Ns and Ng, respectively, we can define 
/S.Nc = Nc — Ns — Ng as an absolute measure of this extra contribu- 
tion.^ Already at first glance, looking at the later part of the simulation, 
we see that this is very significant. In the two lower panels of the Figure, 
we plot histograms of ANc and r = ANc/{Ns + Nc), i.e., the extra 
contribution expressed as a fraction of Ns + Nc- 

The basic observations are as follows. The spiky nature of the grey 
histogram is due to comet showers caused by close stellar encounters 
(we will briefly treat these below). While during the first Gyr the level 
of Ng is generally higher than that of Ns, this situation gets reversed 
after more than two Gyr. Even outside the main showers, Ns is then at 
a somewhat higher level than Ng- The white histogram, showing Nc, 
shares the spikes of the strongest showers, but the contrast between the 
spikes and the background is less than in the grey histogram. Indeed, 
the ANc histogram shows no spikes at all, and the general level does 
not seem to correlate with the stellar injection rate, as illustrated by 
Fig. 3. Therefore, during the later part of the simulation, the r parame- 
ter shows fiuctuations anticorrelated with those of Ns- It reaches a few 
hundred percent, when Ns drops to its lowest levels, but sometimes 
decreases to nearly zero during the peaks of Ns- 

In order to smooth out those fluctuations we present in Table II time 
averages of r over 1 Gyr periods along with the corresponding integrals 
of Nc, Ns and Ng- During the first Gyr the flux of the combined model 



^ Towards the end of our simulation the number of Oort Cloud comets has de- 
creased in all three models but most in the combined one. We then have about 
930 000, 840 000 and 760 000 comets in the tides-only, stars-only and combined 
models, respectively. This means that AA'^o actually underestimates the extra 
contribution of the combined model. 
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Figure 3. ANc versus Ns (plotted on a log scale), where the numbers refer to 
injections of comets into the observable region during intervals of 10 Myr. 

Table II. Number of comets entering the observable region during periods of 
1 Gyr. Model G corresponds to the Galactic tide alone, S to passing stars alone, 
and C to Galactic tide and passing stars together, (r) is the increment from the 
sum of the two first rows (Galactic tide plus passing stars separately) to the 
third row (Galactic tide and passing stars together). 

I Model I [0-1] Gyr [1 - 2] Gyr [2 - 3] Gyr [3 - 4] Gyr [4 - 5] Gyr I 



G 


2128 


797 


481 


307 


248 


S 


1425 


1555 


1030 


717 


511 


c 


3 618 


3141 


2 412 


1733 


1274 


(r) 


1.8% 


33.6% 


59.6% 


69.2% 


67.9% 



is not much larger than the sum of the two fluxes with separate effects, 
and the difference is just a small fraction of the total flux. But toward 
the end the synergy effect of the combined model, as measured by ANq, 
has grown - on the average - to nearly the same level as A''^ -|- Nq- 
During the last Gyr we find that (Nq) is about 2.5 times larger than 
{Ns) in fair agreement with earlier estimates by Heisler et al. (1987) 
and Heisler (1990). After an initial, relatively fast decrease due to the 
emptying of the tidal infeed trajectories, (Nc) continues to decrease 
approximately in proportion to the total number of Oort Cloud comets 
{Noc)i and (Ns) and {No) show similar behaviours. 

Looking in detail at the ANc and r histograms in Fig. 2 for the 
beginning of our simulation, we see that they start from negative values 
and turn into positive ones after ~ 0.5 Gyr. Thus, in the very begin- 
ning, the sum of the separate fluxes is larger than the combined flux. 
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This phenomenon was found by Matese and Lissauer (2002), whose 
calculations were limited to only 5 Myr, and as they explained, it is 
typical of a situation where both tides and stars individually are able 
to fill the loss cone to a high degree. We will discuss this point again 
in Sect. 5. 

The large amount of synergy (r ^ 70%) seen in the later part of 
our simulation is remarkable and indicates that both the tides and the 
stars on their own are seriously inefficient in filling the loss cone. It 
is only by means of the synergy of both effects that we are able to 
explain an important degree of loss cone filling at the current epoch. 
We will look at this closer in Sect. 5 by separating the injection flux 
into different ranges of semi-major axis. For the moment we emphasize 
that treating comet injections from the Oort Cloud in the contemporary 
Solar System simply as a result of the Galactic tide is not a viable idea. 

Already in Sect. 3 we identifled a mechanism that offers a likely 
explanation of the synergy effect, i.e., the repopulation of tidal infeed 
trajectories via stellar encounters. But note in Fig. 2 that the initial 
flux of the model with tides only is not matched by the white areas in 
the later part of the simulation. Thus, even though there is an ongoing 
replenishment of the tidal infeed trajectories due to the randomizing 
effect of stellar encounters, this replenishment is not complete. The 
critical trajectories remain largely depleted, and models that do not take 
this fact into account will overestimate the tidal contribution to the 
injection flux, as well as the efficiency of tides in filling the loss cone. 

The most important synergy mechanism of the Galactic tide and 
stellar perturbations is that the latter are able to repopulate the critical 
phase space trajectories that in the quasi-regular dynamics imposed by 
the tide lead into the loss cone (Dybczyhski, 2002; Fernandez, 2005). 
Our results appear to verify and quantify this picture. But in addition 
we see hints that a different effect is also at work. In Sect. 5 we will 
show that the distribution of 1/a of newly injected, observable comets - 
even during the typical, quiescent periods - has a significant extension 
inside the limit (at a ~ 3 x 10*^ AU), where the tide becomes able 
to feed comets from outside the loss cone into observable orbits. Our 
explanation for this effect is as follows. 

In qualitative terms, when the Galactic tide is in the process of 
injecting a comet into the observable region from the region outside the 
loss cone, and stellar perturbations are added, the latter will sometimes 
aid in decreasing the perihelion distance of the comet (Ag* < 0), and 
sometimes they will counteract the tide (Ag* > 0). To first order, 
the two effects will cancel. But if we consider how much the critical 
value (l/flc) of the inverse semi-major axis for tidal injection into the 
observable region (oc — 3 x 10^ AU) is decreased by a typical positive 
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A(7* or increased by the same typical value of Aq^ in the negative 
direction, we find that the latter effect dominates. Thus, a negative 
Ag* causes a larger gain of comets with a < Uc than the loss of comets 
with a > ttc caused by a positive Ag* of the same size. This holds for 
any nearly uniform distribution of 1/a in the Oort Cloud. 

In mathematical terms, consider the expression for the maximum 
possible decrease of q in one revolution due to the Galactic tide: 

Aq = -Sq^l^z-'/^ (2) 

where z = 1/a and S = 2.8 x 10^"*^^ (Byl, 1986), counting q and a in 
AU. This would hold for a Galactic latitude of perihelion of ±45°. We 
take this favourable orbital orientation as an example, but the following 
arguments apply for any other orientation as well. 

Next, consider a particular value {qp) of the perihelion distance 
preceding the injection into an observable orbit. Using Eq. (2), one 
can write down an approximate condition for the critical value z = zq 
in order to bring the comet into the observable region (g < g^ = 5 AU) : 

qp = qe + Sq^^z^'^'^ (3) 

Eq. (3) defines a unique relation between qp and zq, and by differen- 
tiating one easily finds that qp decreases monotonously with zq, while 
the second derivative is always positive. 

Considering thus an arbitrary point {zo,qp) satisfying Eq. (3), we 
may introduce stellar perturbations by adding a term —Aq^ to the 
right-hand member of Eq. (3), and we can write: 

qp = qi + Sql/^z-'^^-Aq, (4) 

where zi is the new critical value of z. Hence, (zi, gp + Ag*) also satisfies 
Eq. (3). Without the stellar perturbation all comets with q = qp and 
z < Zq are injected into q < qi, and including the stellar perturbation, 
the condition changes into z < zi. From the negative slope of the qpizo) 
relation it is obvious that a negative value of Ag* yields zi-|- > zq, and 
the same positive value yields zi_ < zq. It is also obvious from the 
positive curvature of the graph that zi+ — zq > zq — z\-. If \Aq^\ 
is very small, this difference is negligible, but if it is large enough to 
compete with lAgd, the effect will be important. The latter is indeed 
often the case, when we discuss injections from just outside the loss cone 
{(Ip ~ 15 AU), as has been shown, e.g.^ by Duncan et al. (1987) and 
Fernandez (2005). Therefore, the gain of comets with z > zq occurs over 
a larger interval than the loss with z < zq for a symmetric distribution 
of stellar perturbations. 
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Another issue is the distribution function of z for the Oort Cloud. 
Our simulations start with a probabihty density /(z) oc z^^''^ as ap- 
propriate for an Oort Cloud formed according to the model of Duncan 
et al. (1987). In such a situation there would be more comets per unit 
interval of z at z < zq than at z > zq, and the gain effect would be 
counteracted and possibly turned into a net loss of injected comets. 

However, an interesting result of our simulations is that the gradual 
loss of comets from the Oort Cloud changes the distribution of 1/a. In 
agreement with a recent study by Brasser et al. (2008), we find that the 
loss of comets from the outer parts of the cloud is not compensated by 
diffusion from the inner parts, so that after 4.5 Gyr, when the number 
of comets has decreased from 1 x 10^ to 7.6 x 10^, f{z) has become 
roughly flat over the range from a = 20 000 to 100 000 AU. This shows 
that we have to expect a net gain of newly injected, observable comets 
resulting from a synergy of Aqc and Ag*. Moreover, there should be a 
shift of comets from outside to inside the tidal injection limit - probably 
explaining why we see a significant flux of new comets all the way down 
to a ~ 20 000 AU. 

Although we cannot provide exact numbers, it appears that the 
secondary synergy mechanism due to what we may call "constructive 
interference" of the two effects - even though it certainly exists - is not 
the dominating one. The tentative evidence comes from the relative 
smoothness of the ANc histogram (Fig. 2) and the lack of correlation 
between AA'^c and Ns (Fig. 3). These properties are expected of the 
repopulation of tidal infeed trajectories because of the long response 
time (~ several 10® yr) for tidal infeed on the T^ time scale (Fig. 1). 
But if the constructive interference had been very important, we would 
have expected AA^c' to increase immediately upon an increase of Ns 
- with the caveat that visible peaks of our Ns histogram might arise 
primarily from a temporary infeed of inner Oort Cloud comets with 
z >> Zq, which do not contribute to the interference. Further detailed 
studies are needed to clarify this issue. 

The cometary showers, displayed in Fig. 4 by means of a histogram 
of the injection flux of the combined model with a bin width of only 
10 Myr, are seen to be quite important for the injection of comets from 
the Oort Cloud, as expected and as found by other authors previously 
{e.g., Heisler et al., 1987). We are saving a detailed analysis of those for 
a later paper. At present, we can only remark that the results presented 
here are hard to compare with the treatments of cometary showers by 
Heisler et al. (1987) or Heisler (1990). The first of these papers treated 
only comets with a = 10 000 or 20 000 AU with a full dynamical model 
and then only for a time interval of less than 200 Myr. The second 
gave only a brief account of a simulation for 4.5 Gyr and then only for 
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Figure 4- Same as Fig. 2 (upper panel), but versus time per 10 Myr for the combined 
model only. 



injections from g > 10 to g < 10 AU instead of our requirement that 
the perihelion has to fall substantially from q > 15toQ'<5 AU. 



5. Orbital element distributions of observable comets 

Figure 5 shows the distributions of the opposite of the inverse semi- 
major axis (— 1/a) and the sine of the Galactic latitude of perihelion 
(for clarity we use the absolute value | sin6|) of the comets entering the 
observable region, i.e., heliocentric distance smaller than 5 AU, during 
a typical 170 Myr interval near the end of our simulation, where no 
strong comet showers are registered. We show an average of three such 
periods, i.e., 4.38 - 4.55 Gyr, 4.55 - 4.72 Gyr and 4.80 - 4.97 Gyr. In 
fact, comparing the three data sets, we find a rather good agreement, 
so that tentatively, the expected error of the mean is not very large. 
Three histograms are shown for each quantity: the one in black is for the 
model with Galactic tides only, the grey one is for the model with only 
stellar perturbations, and the white one shows the combined model. 

After more than 4 Gyr the Galactic tides alone are practically only 
able to inject comets into the observable region if a > 50 000 AU, so that 
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Figure 5. Distributions of —1/a, where a is the semi-major axis (top panels) and 
I sin6|, where 6 is the Galactic latitude of perihelion (bottom panels), for the comets 
entering the observable region during 170 Myr near the end of the simulation. When 
present, numbers in the top- left corners of —1/a distribution panels correspond to 
comets with —1/a < — 1 x 10^ AU~ . The left column corresponds to the model 
with Galactic tide alone, the middle column to passing stars alone, and the right 
column to the model with both effects. 



the non-integrable part of the tides may provide new comets into the 
emptied infeed trajectories of the vertical component. Thus the feeble 
flux of new observable comets is strictly confined to the outermost parts 
of the Oort cloud. If only the stellar perturbations are at work, the 
injected comets are almost as few as in the case of the Galactic tides. 
However, the distribution of —1/a shows that the stellar perturbations 
are relatively efficient injectors of comets with semi-major axes in the 
whole range from 25 000 to more than 100 000 AU, and there is some 
marginal infeed all the way into the inner core. Note that this concerns 
a time interval without any strong comet showers. 

When both the processes are at work, the number of comets entering 
the observable zone is 206, about 86% more than the sum of the two 
separate contributions (39 -|- 72). This estimate of r is a bit higher 
than for the entire 1 Gyr interval, listed in Table II, because the three 
intervals have been selected as particularly calm, avoiding even the 
smaller peaks of Ns that can be seen in Fig. 2. We have shown above 
that larger values of Ns lead to smaller values of r. The distribution 
of —1/a is as wide as for the stellar perturbations alone. However, 
the picture has changed, since the additional 86% of the comets are 
strongly concentrated to the interval from —4 x 10~^ to —2 x 10^^ AU~ 
(25 000 < a < 50 000 AU). The local values of AiV^ for the five 1/a 
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intervals (0 - 1), (1 - 2), (2 - 3), (3 - 4) and (4 - 5) x IQ-^ AU^^ 
are —2, —22, +63, +36 and +10, respectively. We will comment on the 
negative values of the first entries in connection with Fig. 7. 

We see that the mechanism of synergy that increases the flux of 
injections in the combined model prefers the range of semi-major axis 
(a > 30 000 AU) where the vertical Galactic tide is able to provide 
the injections, once the relevant trajectories are populated. But there 
is an important extension of the synergy to smaller semi-major axes 
as well, extending at least to a ~ 20 000 AU. We conclude that both 
the above-described synergy mechanisms must be at work. The re- 
population mechanism is obviously important, but the shift to smaller 
semi-major axes can only be explained by the 'constructive interference' 
mechanism. 

Looking at the distributions of jsinfej, indeed the signature of the 
Galactic tide is clearly present in the left diagram and absent in the 
middle one. However, it appears again to some extent in the right- 
hand diagram, where the combined model is presented. Thus we have 
evidence that the synergetic injection of comets in the combined model 
carries an imprint in the latitudes of perihelia similar to that of the 
Galactic tide, though the feature is strongly subdued. In fact, while the 
subdued tidal imprint is consistent with an important role being played 
by the 'constructive interference' synergy mechanism, our combined 
model does not appear to reproduce the observed |sin6| distribution 
of new Oort Cloud comets. An in-depth study of this problem and 
a consideration of ways out of this possible dilemma are deferred to 
future papers. 

The shaping of the b distribution by the Galactic tide was first 
simulated numerically for a realistic Oort Cloud model by Matese and 
Whitman (1989). However, the left panels of Fig. 5 show a behaviour 
that is in stark contrast to their results. Practically all our tidal in- 
jections occur for a > 50 000 AU, where Matese and Whitman (1989) 
found no tidal imprint in the b distribution because of complete loss 
cone filling independent of b. In the light of our results this can be seen 
as an artifact of their assumption of complete randomization of the Oort 
Cloud orbit distribution. Indeed, as we shall find below (Table IV), the 
tidal loss cone filling for 50 000 < a < 100 000 AU towards the end of 
our simulation is far from complete, and therefore we see the imprint 
of the tide in our b distribution. 

In Fig. 6 we show the corresponding distribution of — 1/a and | sin6| 
for the 10 Myr interval from 3.85 to 3.86 Gyr, where Fig. 4 shows that 
the number of observable comets has a high peak due to a moderately 
strong shower. Occurring near the middle of the period, this shower 



spolrev.tex; 16/04/2008; 9:50; p. 20 



21 

dominates the time-integrated injection rate. The trigger is a M5 star 
with an impact parameter d ~ 2 000 AU and a velocity v^, ~ 18 km/s. 
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Figure 6. Same as Fig. 5 but considering the comets that enter the observable region 
during a shower between 3.85 and 3.86 Gyr. The shower is due to a M5 star with 
impact parameter 2 055 AU, velocity 17.7 km/s and mass 0.21 Mq. 

The mid and right-hand —1/a distributions show that, as soon 
as stars are involved, the injection of comets now extends over the 
whole cloud, including an important fraction from the inner core with 
a < 10 000 AU. In fact, the synergy effect is now very strong in 
the range from 10 000 to 20 000 AU, amounting to r > 150%. This 
is unexpected on the basis of both the above-mentioned mechanisms, 
since we are discussing orbits too far inside the tidal injection limit. 
We are instead led to hypothesize a different mechanism. In the present 
case we are comparing the number of comets injected by a particular, 
deeply penetrating star from the mentioned range of semi-major axes 
in the stars-only vs the combined model. In the absence of the Galactic 
tides it is likely that orbits with perihelia close to but outside the loss 
cone have been depleted by the preceding cometary showers, while in 
the combined model the disk tide provides a regular transfer of comets 
into this zone on a Gyr time scale, thus compensating for the losses. 
This means that the synergy works in the opposite sense compared to 
the normal situation outside the showers. During a shower the tides are 
providing the material for injections by the stars, while in the normal 
situation the stars are providing the material for tidal injections. The 
absence of a synergy in the inner core may be explained by the very 
long time scale for tidal torquing, or by a smaller degree of depletion 
of the source region for stellar injections. 
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Naturally, in the stars-only model the shower does not carry any 
signature in the distribution of |sin6|. The combined model does not 
exhibit any significant signature either, but there may nonetheless be a 
slight tendency. In case this is real, it might possibly reveal a somewhat 
more efficient synergy in the 10000 — 20000 AU range for the orbits 
experiencing a faster tidal decrease of q. 
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Figure 7. Same as Fig. 5 but considering the comets entering the observable region 
during 170 Myr near the beginning of the simulation. 

Let US now consider the situation at the beginning of the simulation, 
before the tides have had the time to completely empty the tidal infeed 
trajectories in the outer part of the cloud. The results are shown in 
Fig. 7 for a period between 0.63 and 0.80 Gyr, when no strong showers 
are noted. The number of comets entering into the observable region 
is 282, 128 and 463 for the tides, the stars, and the combined model, 
respectively. The action of the tides is still quite strong, since the infeed 
trajectories in the interval 20 000 < a < 50 000 AU are not yet severely 
depleted (cf. Fig. 1). Therefore the net synergy effect amounts to only 
13%. 

The local ^Nc values for the same five l/a intervals as we discussed 
in connection with Fig. 5 are in this case —5, —66, +35, +61 and 
+ 19. We do not see any significant synergy effect for more tightly 
bound orbits. The distribution of positive synergy over the 2 — 5 x 
10^^ AU^^ range is similar to that of Fig. 5, and our conclusions about 
the relevance of the two mechanisms are the same. Note that in both 
cases we see negative A.Nc values in the two outermost l/a ranges 
(0 — 2 X 10~^ AU^^). The fundamental reason is the one discussed 
above in order to explain the negative r values in the very beginning of 
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the simulation, i.e., a saturation effect of the loss cone when both injec- 
tion effects individually are able to cause a large degree of filling. The 
distributions of |sin6| exhibit the same features as in Fig. 5 and lead 
to the same conclusion: when both tides and stars act in synergy, the 
signature of the Galactic tide may be seen but appears quite marginal. 
Since we noted in Sect. 4 that the replenishment of the tidal infeed 
trajectories by stellar perturbations is not complete during the later 
part of our simulation, an obvious consequence is that the filling of the 
loss cone for the relevant semi-major axes cannot be complete either. To 
quantify this statement, we consider the rate of perihelion passages h 
with q < qo as a function of the semi-major axis, assuming complete loss 
cone filling and a completely thermalized Oort cloud (Hills, 1981; Bailey 
and Stagg, 1988): 



h{a) = Noc ■ f{a) a 

a 



-3/2 



(5) 



Here, iVoc is the number of comets in the entire Oort cloud (ini- 
tially 10^ in our model) and f{a) the frequency function describing the 
distribution of semi-major axes: /(a) oc a~^'^ initially in our model. 

We have computed injection rates in the combined model for the 
three time intervals of Figs. 5-7 using Eq. (5) and finding the inte- 
grals / Nocf{0')0'~^ da over different ranges of — 1/a directly from the 
simulation output at neighbouring moments. The calculation of these 
integrals is done by simply adding the values of a~^ " of all the comets 
found in the relevant range. Multiplying by the length of each interval 
and putting qq = 5 AU, we find the numbers iVcomp listed in Table HI. 

Table III. Numbers of comet injections during the time intervals of Figs. 5-7 for different 
ranges of inverse semi-major axis, as computed from Eq. (5) and found from our simu- 
lation of the combined model. The ratio of simulated to computed number, expressed in 
percent, is also listed in each case. 

A(l/a) Beginning End, sho\vcrs End, quiescent 

(10^^ AU"^) ATcomp iv~ 



/ic 



Nsin 



/lo 



Afein 



(0-2) 


106 


111 


~ 100% 


3.2 


2 


~ 60% 


48 


40 


83% 


(2-3) 


303 


189 


62% 


10 


6 


~ 60% 


160 


83 


52% 


(3-4) 


616 


111 


18% 


24 


2 


~ 8% 


367 


46 


13% 


(4 - 5) 


1044 


23 


2.2% 


44 


6 


14% 


692 


14 


2.0% 


(5 - 10) 


15600 


23 


0.15% 


740 


29 


3.9% 


12100 


13 


0.11% 


> 10 


672000 


6 


0.0009% 


37200 


31 


0.08% 


626000 


10 


0.0016% 



The numbers of comet injections for each 1/a range and each time 
interval can be read off from the figures, and they are listed as Nsim in 
the Table along with the ratios A'sim/A''comp; which give the filling factor 
of the observable part of the loss cone (/ic)- We find that this factor is 
close to 100% in the beginning of the simulation for a > 50 000 AU and 
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Table IV. Filling factors for the observable part of the loss cone, com- 
puted for different ranges of semi-major axis and separately for the three 
dynamical models (tides-only, stars-only, and combined). 

I A(l/a) I 



/ic (Beginning) 



Stellar Combined Tidal 



/ic (End) 



Stella 



(0-1) 


~ 400% 


~ 100% 


~ 200% 


~ 100% 


~ 100% 


~ 100% 


(1-2) 


86% 


60% 


~ 100% 


30% 


45% 


78% 


(2-3) 


36% 


10% 


62% 


0.6% 


10% 


52% 


(a -4) 


6.5% 


1.1% 


18% 


- 


2.3% 


13% 


(4-5) 


0.09% 


0.3% 


2.2% 


- 


0.5% 


2.0% 


(5 - 10) 


- 


0.1% 


0.15% 


- 


0.06% 


0.11% 


> 10 


- 


0.0006% 


0.0009% 


- 


0.0008% 


0.0016% 



remains > 80% for such semi-major axes even towards the end during 
quiescent periods. But the factor drops rapidly with decreasing a to 
values near 2% at a ~ 20 000 AU. These results may be compared with 
those of Heisler (1990), who used a similar procedure for deriving fi^. 
She did not consider semi-major axes a > 40 000 AU, and inside this 
limit we find somewhat smaller filling factors than she did, consistent 
with the fact that we use a lower value for the Galactic mid-plane 
density and somewhat higher stellar velocities. 

Note that the filling factors have decreased somewhat, when we 
compare the final quiescent period with the initial one. Except in the 
outermost parts of the Oort Cloud, there is always a depletion of comets 
in the regions of phase space near the tidal infeed trajectories and in 
the vicinity of the loss cone, and this depletion grows slowly with time. 

The numbers A'sim found for the shower period are too small to 
be statistically useful for the outer parts of the cloud, and the filling 
factors listed are very uncertain. However, we see an obvious effect in 
the inner parts, when comparing /ic with the corresponding values of 
quiescent periods. The shower increases fic by factors ~ 20 — 100, and 
thus the overall flux exhibits the peak seen in Fig. 4 due to comets with 
a < 20 000 AU. 

We have already made the remark that neither ANc nor r provides 
a fully satisfactory measure of the synergy effect, because they do not 
account for the difference of the number of Oort Cloud comets between 
different dynamical models ~ especially towards the end of the simula- 
tion. After 5 Gyr the total number of comets in the combined model is 
only ~ 80% of that in the tides-only model, and if we concentrate on 
comets with 50 000 < a < 100 000 AU where the losses are the largest, 
the ratio of the two models is only 35%. In order to compensate for 
such effects we have computed the fic parameter separately for the 
three models and for all the ranges of 1/a, and we present the results 
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in Table IV. The time periods referred to are the quiescent periods of 
Figs. 5 and 7. 

The outermost energy range is empty in all models, when the sim- 
ulation starts, but it gets populated quickly - at least when stars are 
involved. We interpret the very large value of /ic in the tides-only model 
at the beginning as evidence that the radial tide has not extracted 
comets into this energy range in a uniform manner, so that our as- 
sumption of thermalization when deriving A/'comp is not justified. To a 
lesser extent this appears to be true also in the combined model, where 
stars have extracted many more comets. It is likely that this extraction 
too - at the early time in question - has not populated all the angular 
momenta in a thermalized fashion. However, the statistics is too poor 
to be confident about such conclusions. In any case, the loss cone filling 
is extremely efficient for all models, thus explaining the negative values 
of ANc. 

For the next energy range we see the saturation effect again, espe- 
cially in the beginning. At the end, the value of fic in the combined 
model is close to the sum of those in the other two models. Hence there 
is no apparent synergy in this case, but probably a real synergy has 
been concealed by the saturation effect. In any case the large negative 
value of ANc results entirely from the small number of comets in the 
combined model, as discussed above. In the next three energy ranges 
(20 000 < a < 50 000 AU) we see that /ic in the combined model is much 
larger than the sum of the two other entries, and for a < 20 000 AU the 
effect continues: adding the tides to the stars increases the loss cone 
filling by a factor 1.5 — 2. 



6. Discussion and conclusions 

We have simulated the evolution of the Oort cloud over 5 Gyr, using for 
initial conditions a relaxed model with a distribution of semi-major axis 
/(a) oc a-i-^ within the interval 3 000-100 000 AU. This model is based 
on the results of simulation of Oort Cloud formation and evolution by 
Duncan et al. (1987). However, we do not find this to be a steady 
distribution. More comets are lost from the outer parts of the cloud 
than can be replaced from inside, so that our model cloud evolves into 
a distribution close to f{a) oc a~^ ~ i.e., flat in 1/a. 

Our dynamical model has two main limitations. We do not treat 
encounters with very massive Galactic perturbers, such as star clusters 
or Giant Molecular Cloud complexes, the justiflcation being that they 
occur so rarely that the current Solar System is unlikely to feel the 
direct reverberations of any such encounter, and that even if they 
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modify the structure of the Oort cloud, our interest is not primarily 
in its dynamical history but rather in the way stars and Galactic tides 
currently interact when injecting observable comets. 

Moreover, we do not treat planetary perturbations in any direct 
manner. Like most previous investigators {e.g., Heisler, 1990) we use a 
loss cone defined by a limiting perihelion distance (in our case, 15 AU) 
outside which no planetary effects are included and inside which all 
comets are considered lost from the cloud through perturbations by 
Jupiter and Saturn. In terms of "transparency" of the planetary system 
(Dybczyhski and Pr§tka, 1997; Dybczyhski, 2005), our model is com- 
pletely opaque (P = 1). This means that we are limiting our attention 
to a subset of the observed population of "new" Oort cloud comets, i.e., 
those that have jumped directly from q > 15 AU into their observed 
orbits with q < 5 AU. We are neglecting the rest of the population, 
which consists of comets that passed perihelia in the outer part of the 
loss cone without being perturbed away before arriving into observable 
orbits. We are also neglecting a possible contribution to the observed 
new comets by a "leakage" from the scattered disk (Levison et al., 
2006). Therefore we prefer not to draw any conclusions in this paper 
regarding the total number of Oort Cloud comets or the exact values 
of the filling factors. Nor do we claim to make any prediction on the 
detailed shape of the 1/a distribution of new Oort Cloud comets, until 
we have included the planetary perturbations in a realistic manner. 

We have shown that the concept of tidal and stellar torquing time 
scales (Duncan et al., 1987) gives a very incomplete picture of the 
speed of comet injection, whether it may concern Galactic tides or 
stellar encounters. The distribution of injection times is largely shaped 
by other effects - like comet showers or the repopulation of the emptied 
infeed trajectories of the disk tide due to the non-integrable part of the 
tides or stellar perturbations. 

We have also shown how - for semi-major axes large enough for 
the tide to populate observable orbits - the regions of the phase space 
occupied by trajectories leading into the loss cone get depleted during 
the first Gyr of Oort cloud evolution. This would leave little chance for 
the tide to produce a significant number of observable comets at the 
current time, were it not for the capability of stellar perturbations to 
replenish the tidal infeed trajectories. 

We have indeed demonstrated that, during the later parts of our 
simulation, there is a very important synergy effect of the Galactic tide 
and stellar perturbations such that the combined injection rate is on 
the average ~ 70% larger than that of the stars alone plus that of 
the tide alone. This synergy is strongest for semi-major axes between 
~ 20 000 — 50 000 AU but continues all the way into the inner core. 
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During comet showers the synergy effect in the outer parts of the cloud 
practically disappears, but the one affecting the inner parts becomes 
very important. 

We have identified two mechanisms for the synergy during quies- 
cent periods in the outer parts of the Oort Cloud. One is that the 
stellar perturbations provide a supply of new comets that replenishes 
the depleted tidal infeed trajectories, and the other is that the gain 
of comet injections, when stellar perturbations decrease the perihelion 
distance, dominates over the loss caused by opposing perturbations. For 
the synergy of the inner cloud we hypothesize that the Galactic tides 
provide the material for stellar injections by slowly feeding the region 
of phase space in the vicinity of the loss cone. Thus, the general picture 
spawned by our results is that injection of comets from the Oort Cloud 
is essentially to be seen as a team work involving both tides and stars. 
It appears meaningless to rank the two effects in terms of strength or 
efficiency. 

Indeed, for the smaller semi-major axes the Galactic tide does not 
dominate the injection of comets, contrary to the conclusions of Heisler 
et al. (1987) and Heisler (1990).^ It only contributes to a synergy with 
stellar perturbations, and without the stars one would not have any 
injections of comets with o ^ 20 000 AU. 

The distribution of Galactic latitudes of perihelia of the observable 
comets exhibits a maximum for | sin6[ ~ 0.5 as expected in the tides- 
only model, but in the combined model this feature can hardly be seen 
at all. The tides form part of the synergetic injection, but their imprint 
is largely washed out by the stellar contribution. But, likely due to 
the role of the tides in helping the stars to create comet showers, the 
pattern can be seen at least as clearly during a shower as during the 
quiescent periods. Therefore, it tentatively appears that the shape of 
the observed b distribution can not be used to indicate whether we 
are experiencing any shower at present. However, since none of our 
model distributions appears to agree with the observed one, we have to 
defer any conclusions to future papers. It may be interesting to see, for 
instance, if the leakage from the scattered disk into the Oort Cloud with 
an ensuing direct transfer into observable new comets may alleviate the 
problem. 



^ The main reason for this discrepancy is that the Heisler papers considered 
injections into the loss cone - mainly by slight perturbations of q across the limiting 
value Qc — 10 AU - while we consider large jumps from q > 15 AU into the 
observable region with q < 5 AU. Interestingly, Heisler et al. (1991) commented that 
the injection into orbits with a <, 20 000 AU and g < 2 AU is indeed dominated 
by stellar perturbations. 
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We have measured the filhng of the observable part of the loss cone 
by comparing our simulated injection rates for different intervals of 
semi-major axis with the rates of observable perihelion passages {q < 
5 AU) computed for a completely thermalized distribution of cometary 
orbits involving a filled loss cone. The deficiency of our simulated rate 
likely reflects not only a lack of comets in the loss cone but a general 
depletion in a wider phase space region in its vicinity, as remarked by 
Heisler (1990). Our results can be compared with hers, and in contrast 
to her inference that /ic may level out at ~ 60% for a ^ 30 000 AU, 
we find an average filling factor during quiescent periods in the current 
Solar System, which drops steadily from ~ 100% at a > 100 000 AU 
to 1% or less at a < 20 000 AU in the combined model. However, there 
are important differences between the two investigations, one being that 
she simulated a much shorter time period than we do, and in addition 
our parameters for the Galactic tides and stellar encounters also differ 
from hers. 

In agreement with Weissman and Hut (1986), we find that cometary 
showers do not significantly increase the loss cone filling at large semi- 
major axes. However, near 25 000 AU there is an abrupt change into the 
regime of the inner cloud, where the filling factor increases by orders 
of magnitude during moderate to strong events. The showers of course 
involve direct injections by the passing stars, but the synergy with the 
Galactic tide is as important as during quiescent periods. 
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